3 Equivalent Water Thickness/Canopy Water Content from Imaging Spectroscopy Data

Summary

In this notebook we will explore how Equivalent Water Thickness (EWT) or Canopy Water Content (CWC) can be calculated from the Earth Surface Mineral Dust Source Investigation (EMIT) L2A Reflectance Product, then we will apply this knowledge to calculate CWC over the Jack and Laura Dangermond Preserve located near Santa Barbara, CA.

Background

Equivalent Water Thickness (EWT) is the predicted thickness or absorption path length in centimeters (cm) of water that would be required to yield an observed spectra. In the context of vegetation, this is equivalent to canopy water content (CWC) in g/cm^2 because a cm^3 of water has a mass of 1g.

CWC can be derived from surface reflectance spectra because they provide information about the composition of the target, including water content. Reflectance is the fraction of incoming solar radiation reflected by Earth’s surface. Different materials reflect varying proportions of radiation based upon their chemical composition and physical properties, giving materials their own unique spectral signature or fingerprint. In particular, liquid water causes characteristic absorption features to appear in the near-infrared wavelengths of the solar spectrum, which enables an estimation of its content.

CWC correlates with vegetation type and health, as well as wildfire risk. The methods used here to calculate CWC are based on the ISOFIT python package. The Beer-Lambert physical model used to calculate CWC is described in Green et al. (2006) and Bohn et al. (2020). It uses wavelength-dependent absorption coefficients of liquid water to determine the absorption path length as a function of absorption feature depth. Of note, this model does not account for multiple scattering effects within the canopy and may result in oversestimation of CWC (Bohn et al., 2020).

The Jack and Laura Dangermond Preserve and its surrounding lands are one one of the last “wild coastal” regions in Southern California. The preserve is over 24,000 acres and consists of several ecosystem types and is home to over 600 plant species and over 200 wildlife species.

More about the EMIT mission and EMIT products.

References

• Shrestha, Rupesh. 2023. Equivalent water thickness/canopy water content from hyperspectral data. Jupyter Notebook. Oak Ridge National Laboratory Distributed Active Archive Center. https://github.com/rupesh2/ewt_cwc/tree/main
• Bohn, N., L. Guanter, T. Kuester, R. Preusker, and K. Segl. 2020. Coupled retrieval of the three phases of water from spaceborne imaging spectroscopy measurements. Remote Sensing of Environment 242:111708. https://doi.org/10.1016/j.rse.2020.111708
• Green, R.O., T.H. Painter, D.A. Roberts, and J. Dozier. 2006. Measuring the expressed abundance of the three phases of water with an imaging spectrometer over melting snow. Water Resources Research 42:W10402. https://doi.org/10.1029/2005WR004509
• Thompson, D.R., V. Natraj, R.O. Green, M.C. Helmlinger, B.-C. Gao, and M.L. Eastwood. 2018. Optimal estimation for imaging spectrometer atmospheric correction. Remote Sensing of Environment 216:355–373. https://doi.org/10.1016/j.rse.2018.07.003

Requirements - NASA Earthdata Account
- No Python setup requirements if connected to the workshop cloud instance!
- Local Only - Set up Python Environment. See setup_instructions.md in the /setup/ folder
- Local Only - Downloaded necessary files. This is done at the end of the 01_Finding_Concurrent_Data notebook.

Learning Objectives
- Calculate CWC of a single pixel
- Calculate CWC of an ROI

Tutorial Outline

3.1 Setup
3.2 Opening EMIT Data
3.3 Extracting Reflectance of a Pixel
3.4 Calculating CWC
3.4.1 Single Point
3.4.2 DataFrame of Points
3.5 Applying Inversion in Parallel Across an ROI

# Install ray on 2i2c
!pip install "ray[default]"

Requirement already satisfied: ray[default] in /srv/conda/envs/notebook/lib/python3.10/site-packages (2.8.1)
Requirement already satisfied: click>=7.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (8.1.7)
Requirement already satisfied: filelock in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (3.12.4)
Requirement already satisfied: jsonschema in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (4.19.1)
Requirement already satisfied: msgpack<2.0.0,>=1.0.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.0.5)
Requirement already satisfied: packaging in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (23.1)
Requirement already satisfied: protobuf!=3.19.5,>=3.15.3 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (4.25.1)
Requirement already satisfied: pyyaml in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (6.0.1)
Requirement already satisfied: aiosignal in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.3.1)
Requirement already satisfied: frozenlist in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.4.0)
Requirement already satisfied: requests in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (2.31.0)
Requirement already satisfied: numpy>=1.19.3 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.24.4)
Requirement already satisfied: aiohttp>=3.7 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (3.8.5)
Requirement already satisfied: aiohttp-cors in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (0.7.0)
Requirement already satisfied: colorful in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (0.5.5)
Requirement already satisfied: py-spy>=0.2.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (0.3.14)
Requirement already satisfied: gpustat>=1.0.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.1.1)
Requirement already satisfied: opencensus in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (0.11.3)
Requirement already satisfied: pydantic<2 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.10.13)
Requirement already satisfied: prometheus-client>=0.7.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (0.17.1)
Requirement already satisfied: smart-open in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (6.4.0)
Requirement already satisfied: virtualenv<20.21.1,>=20.0.24 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (20.21.0)
Requirement already satisfied: grpcio>=1.42.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from ray[default]) (1.59.3)
Requirement already satisfied: attrs>=17.3.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from aiohttp>=3.7->ray[default]) (23.1.0)
Requirement already satisfied: charset-normalizer<4.0,>=2.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from aiohttp>=3.7->ray[default]) (3.2.0)
Requirement already satisfied: multidict<7.0,>=4.5 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from aiohttp>=3.7->ray[default]) (6.0.4)
Requirement already satisfied: async-timeout<5.0,>=4.0.0a3 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from aiohttp>=3.7->ray[default]) (4.0.3)
Requirement already satisfied: yarl<2.0,>=1.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from aiohttp>=3.7->ray[default]) (1.9.2)
Requirement already satisfied: nvidia-ml-py>=11.450.129 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from gpustat>=1.0.0->ray[default]) (12.535.133)
Requirement already satisfied: psutil>=5.6.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from gpustat>=1.0.0->ray[default]) (5.9.5)
Requirement already satisfied: blessed>=1.17.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from gpustat>=1.0.0->ray[default]) (1.20.0)
Requirement already satisfied: typing-extensions>=4.2.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pydantic<2->ray[default]) (4.8.0)
Requirement already satisfied: distlib<1,>=0.3.6 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from virtualenv<20.21.1,>=20.0.24->ray[default]) (0.3.7)
Requirement already satisfied: platformdirs<4,>=2.4 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from virtualenv<20.21.1,>=20.0.24->ray[default]) (3.10.0)
Requirement already satisfied: jsonschema-specifications>=2023.03.6 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from jsonschema->ray[default]) (2023.7.1)
Requirement already satisfied: referencing>=0.28.4 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from jsonschema->ray[default]) (0.30.2)
Requirement already satisfied: rpds-py>=0.7.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from jsonschema->ray[default]) (0.10.3)
Requirement already satisfied: opencensus-context>=0.1.3 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from opencensus->ray[default]) (0.1.3)
Requirement already satisfied: google-api-core<3.0.0,>=1.0.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from opencensus->ray[default]) (2.14.0)
Requirement already satisfied: idna<4,>=2.5 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from requests->ray[default]) (3.4)
Requirement already satisfied: urllib3<3,>=1.21.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from requests->ray[default]) (1.26.15)
Requirement already satisfied: certifi>=2017.4.17 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from requests->ray[default]) (2023.7.22)
Requirement already satisfied: wcwidth>=0.1.4 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from blessed>=1.17.1->gpustat>=1.0.0->ray[default]) (0.2.6)
Requirement already satisfied: six>=1.9.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from blessed>=1.17.1->gpustat>=1.0.0->ray[default]) (1.16.0)
Requirement already satisfied: googleapis-common-protos<2.0.dev0,>=1.56.2 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (1.61.0)
Requirement already satisfied: google-auth<3.0.dev0,>=2.14.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (2.25.1)
Requirement already satisfied: cachetools<6.0,>=2.0.0 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from google-auth<3.0.dev0,>=2.14.1->google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (5.3.2)
Requirement already satisfied: pyasn1-modules>=0.2.1 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from google-auth<3.0.dev0,>=2.14.1->google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (0.3.0)
Requirement already satisfied: rsa<5,>=3.1.4 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from google-auth<3.0.dev0,>=2.14.1->google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (4.9)
Requirement already satisfied: pyasn1<0.6.0,>=0.4.6 in /srv/conda/envs/notebook/lib/python3.10/site-packages (from pyasn1-modules>=0.2.1->google-auth<3.0.dev0,>=2.14.1->google-api-core<3.0.0,>=1.0.0->opencensus->ray[default]) (0.5.1)

3.1 Setup

Import the required Python libraries.

# Import Packages
import os
import glob
import earthaccess
import math
import numpy as np
import xarray as xr
from osgeo import gdal
import rasterio as rio
import rioxarray as rxr
from matplotlib import pyplot as plt
import hvplot.xarray
import hvplot.pandas
import holoviews as hv
import pandas as pd
import geopandas as gp
import sys
from modules.emit_tools import emit_xarray, ortho_xr
from modules.ewt_calc import calc_ewt
from scipy.optimize import least_squares

3.2 Open an EMIT File

EMIT L2A Reflectance Data are distributed in a non-orthocorrected spatially raw NetCDF4 (.nc) format consisting of the data and its associated metadata. To work with this data, we will use the emit_xarray function from the emit_tools.py module included in the repository.

Set a filepath for an EMIT L2A reflectance file and open using the emit_xarray function.

fp = '../../shared/2023-VITALS-Workshop-AGU/data/EMIT_L2A_RFL_001_20230401T203751_2309114_002.nc'

ds = emit_xarray(fp,ortho=True)

Similarly to what we did in the previous notebook, we can plot a single band to get an idea of where our scene is. This is the same scene we processed earlier.

emit_layer = ds.sel(wavelengths=850,method='nearest')
emit_layer.hvplot.image(cmap='viridis',geo=True, tiles='ESRI', frame_width=720,frame_height=405, alpha=0.7, fontscale=2).opts(
    title=f"{emit_layer.wavelengths:.3f} {emit_layer.wavelengths.units}", xlabel='Longitude',ylabel='Latitude')

3.3 Extracting Reflectance of a Single Pixel

We can mask out the -.01 values used to represent the region of the spectra with strong atmospheric water vapor absorption features.

ds['reflectance'].data[:,:,ds['good_wavelengths'].data==0] = np.nan

Retrieve the spectra from a sample point by providing a latitude and longitude along with a method using the sel function. This will select the pixel closest to the provided coordinates.

point = ds.sel(latitude=34.5399,longitude=-120.3529, method='nearest')
point

<xarray.Dataset>
Dimensions:           (wavelengths: 285)
Coordinates:
  * wavelengths       (wavelengths) float32 381.0 388.4 ... 2.486e+03 2.493e+03
    fwhm              (wavelengths) float32 ...
    good_wavelengths  (wavelengths) float32 1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0
    latitude          float64 34.54
    longitude         float64 -120.4
    elev              float32 130.3
    spatial_ref       int64 0
Data variables:
    reflectance       (wavelengths) float32 0.01718 0.01754 ... 0.02635 0.02455
Attributes: (12/40)
    ncei_template_version:             NCEI_NetCDF_Swath_Template_v2.0
    summary:                           The Earth Surface Mineral Dust Source ...
    keywords:                          Imaging Spectroscopy, minerals, EMIT, ...
    Conventions:                       CF-1.63
    sensor:                            EMIT (Earth Surface Mineral Dust Sourc...
    instrument:                        EMIT
    ...                                ...
    spatial_ref:                       GEOGCS["WGS 84",DATUM["WGS_1984",SPHER...
    geotransform:                      [-1.20992414e+02  5.42232520e-04 -0.00...
    day_night_flag:                    Day
    title:                             EMIT L2A Estimated Surface Reflectance...
    granule_id:                        EMIT_L2A_RFL_001_20230401T203751_23091...
    Orthorectified:                    True

xarray.Dataset

• Dimensions:
  • wavelengths: 285
• Coordinates: (7)
  • wavelengths
    (wavelengths)
    float32
    381.0 388.4 ... 2.486e+03 2.493e+03

    long_name :

    Wavelength Centers

    units :

    nm

    array([ 381.00558,  388.4092 ,  395.81583, ..., 2478.153  , 2485.5386 ,
           2492.9238 ], dtype=float32)

  • fwhm
    (wavelengths)
    float32
    ...

    long_name :

    Full Width at Half Max

    units :

    nm

    [285 values with dtype=float32]

  • good_wavelengths
    (wavelengths)
    float32
    1.0 1.0 1.0 1.0 ... 1.0 1.0 1.0 1.0

    long_name :

    Wavelengths where reflectance is useable: 1 = good data, 0 = bad data

    units :

    unitless

    array([1., 1., 1., ..., 1., 1., 1.], dtype=float32)

  • latitude
    ()
    float64
    34.54

    long_name :

    Latitude (WGS-84)

    units :

    degrees north

    array(34.5401589)

  • longitude
    ()
    float64
    -120.4

    long_name :

    Longitude (WGS-84)

    units :

    degrees east

    array(-120.35285082)

  • elev
    ()
    float32
    130.3

    long_name :

    Surface Elevation

    units :

    m

    array(130.26474, dtype=float32)

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    array(0)

• Data variables: (1)
  • reflectance
    (wavelengths)
    float32
    0.01718 0.01754 ... 0.02635 0.02455

    long_name :

    Surface Reflectance

    units :

    unitless

    array([0.01718401, 0.01754001, 0.01789793, 0.01825578, 0.01887722,
           0.01977932, 0.02089375, 0.02213017, 0.02327097, 0.02443924,
           0.02553891, 0.02647934, 0.02731651, 0.02806388, 0.03012716,
           0.0303051 , 0.03147712, 0.03389502, 0.03716784, 0.04125404,
           0.04527088, 0.04904396, 0.0513751 , 0.05365279, 0.05416293,
           0.05240294, 0.05057795, 0.04905905, 0.04850443, 0.04918665,
           0.04959129, 0.0485189 , 0.04756864, 0.04709126, 0.04778447,
           0.04722033, 0.04627838, 0.04591437, 0.04385856, 0.04450385,
           0.04487438, 0.04721367, 0.05855059, 0.07728622, 0.09845995,
           0.11936122, 0.14146721, 0.16280495, 0.17837006, 0.1891758 ,
           0.19495177, 0.19879304, 0.20209542, 0.20503408, 0.20749208,
           0.20980741, 0.21209788, 0.21429807, 0.2161592 , 0.21820003,
           0.2202218 , 0.2222354 , 0.22431302, 0.22625907, 0.22809783,
           0.2298785 , 0.23154566, 0.23314181, 0.23455934, 0.2356244 ,
           0.23671815, 0.2378939 , 0.2389565 , 0.23990753, 0.23933375,
           0.238821  , 0.23855045, 0.23721987, 0.23602907, 0.23627467,
           0.23715375, 0.2386118 , 0.2409895 , 0.24408886, 0.24704838,
           0.25001273, 0.25301102, 0.25599742, 0.258825  , 0.26146266,
           0.26378876, 0.2659258 , 0.2676973 , 0.26931012, 0.27055526,
           0.27149647, 0.2718703 , 0.27186936, 0.2714865 , 0.2693018 ,
    ...
                  nan,        nan,        nan,        nan,        nan,
                  nan,        nan,        nan,        nan,        nan,
                  nan,        nan,        nan,        nan,        nan,
                  nan,        nan,        nan,        nan,        nan,
                  nan,        nan,        nan, 0.05113549, 0.05279826,
           0.05452276, 0.05609974, 0.05744791, 0.05934713, 0.06077929,
           0.06164304, 0.06214674, 0.06242346, 0.06241743, 0.06247758,
           0.06254932, 0.062612  , 0.06495836, 0.06524801, 0.06746189,
           0.06759153, 0.06831566, 0.06865513, 0.06930527, 0.07104894,
           0.07243683, 0.07320293, 0.07459541, 0.07555467, 0.07578479,
           0.07651836, 0.07675724, 0.07744075, 0.07785628, 0.0794607 ,
           0.07826924, 0.0792022 , 0.07868779, 0.07894044, 0.07777485,
           0.07542186, 0.07294796, 0.07028439, 0.06905952, 0.06607912,
           0.06531868, 0.06337301, 0.0620889 , 0.05987899, 0.05752026,
           0.05735842, 0.05715922, 0.05624693, 0.05384592, 0.05267731,
           0.05235471, 0.05194782, 0.05061792, 0.05089531, 0.05007195,
           0.04892436, 0.04694235, 0.04603988, 0.0461389 , 0.04569674,
           0.0415682 , 0.04030755, 0.03803281, 0.03684125, 0.03320232,
           0.03311669, 0.03205504, 0.02982526, 0.0263453 , 0.02454741],
          dtype=float32)

• Indexes: (1)
  • wavelengths
    PandasIndex

    PandasIndex(Float64Index([ 381.0055847167969,  388.4092102050781,  395.8158264160156,
                  403.22540283203125, 410.63800048828125,  418.0535888671875,
                   425.4721374511719,  432.8927001953125,  440.3172607421875,
                   447.7427978515625,
                  ...
                   2426.440185546875,  2433.830322265625,   2441.21826171875,
                     2448.6064453125,  2455.994384765625,  2463.381591796875,
                   2470.767822265625,  2478.153076171875,   2485.53857421875,
                      2492.923828125],
                 dtype='float64', name='wavelengths', length=285))

• Attributes: (40)

  ncei_template_version :

  NCEI_NetCDF_Swath_Template_v2.0

  summary :

  The Earth Surface Mineral Dust Source Investigation (EMIT) is an Earth Ventures-Instrument (EVI-4) Mission that maps the surface mineralogy of arid dust source regions via imaging spectroscopy in the visible and short-wave infrared (VSWIR). Installed on the International Space Station (ISS), the EMIT instrument is a Dyson imaging spectrometer that uses contiguous spectroscopic measurements from 410 to 2450 nm to resolve absoprtion features of iron oxides, clays, sulfates, carbonates, and other dust-forming minerals. During its one-year mission, EMIT will observe the sunlit Earth's dust source regions that occur within +/-52° latitude and produce maps of the source regions that can be used to improve forecasts of the role of mineral dust in the radiative forcing (warming or cooling) of the atmosphere.\n\nThis file contains L2A estimated surface reflectances and geolocation data. Reflectance estimates are created using an Optimal Estimation technique - see ATBD for details. Reflectance values are reported as fractions (relative to 1). Geolocation data (latitude, longitude, height) and a lookup table to project the data are also included.

  keywords :

  Imaging Spectroscopy, minerals, EMIT, dust, radiative forcing

  Conventions :

  CF-1.63

  sensor :

  EMIT (Earth Surface Mineral Dust Source Investigation)

  instrument :

  EMIT

  platform :

  ISS

  institution :

  NASA Jet Propulsion Laboratory/California Institute of Technology

  license :

  https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/

  naming_authority :

  LPDAAC

  date_created :

  2023-04-06T02:26:24Z

  keywords_vocabulary :

  NASA Global Change Master Directory (GCMD) Science Keywords

  stdname_vocabulary :

  NetCDF Climate and Forecast (CF) Metadata Convention

  creator_name :

  Jet Propulsion Laboratory/California Institute of Technology

  creator_url :

  https://earth.jpl.nasa.gov/emit/

  project :

  Earth Surface Mineral Dust Source Investigation

  project_url :

  https://emit.jpl.nasa.gov/

  publisher_name :

  NASA LPDAAC

  publisher_url :

  https://lpdaac.usgs.gov

  publisher_email :

  lpdaac@usgs.gov

  identifier_product_doi_authority :

  https://doi.org

  flight_line :

  emit20230401t203751_o09114_s000

  time_coverage_start :

  2023-04-01T20:37:51+0000

  time_coverage_end :

  2023-04-01T20:38:03+0000

  software_build_version :

  010610

  software_delivery_version :

  010610

  product_version :

  V001

  history :

  PGE Run Command: {python /beegfs/store/emit/ops/repos/emit-sds-l2a/spectrum_quality.py /tmp/emit/ops/emit20230401t203751_emit.L2AReflectance_20230405t233102/output/emit20230401t203751_rfl /tmp/emit/ops/emit20230401t203751_emit.L2AReflectance_20230405t233102/output/emit20230401t203751_rfl_quality.txt}, PGE Input Files: {radiance_file=/beegfs/store/emit/ops/data/acquisitions/20230401/emit20230401t203751/l1b/emit20230401t203751_o09114_s000_l1b_rdn_b0106_v01.img, pixel_locations_file=/beegfs/store/emit/ops/data/acquisitions/20230401/emit20230401t203751/l1b/emit20230401t203751_o09114_s000_l1b_loc_b0106_v01.img, observation_parameters_file=/beegfs/store/emit/ops/data/acquisitions/20230401/emit20230401t203751/l1b/emit20230401t203751_o09114_s000_l1b_obs_b0106_v01.img, surface_model_config=/beegfs/store/emit/ops/repos/emit-sds-l2a/surface/surface_20221020.json}

  crosstrack_orientation :

  as seen on ground

  easternmost_longitude :

  -119.71545648976226

  northernmost_latitude :

  34.7990749308955

  westernmost_longitude :

  -120.992414074966

  southernmost_latitude :

  33.77100207248943

  spatialResolution :

  0.000542232520256367

  spatial_ref :

  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

  geotransform :

  [-1.20992414e+02 5.42232520e-04 -0.00000000e+00 3.47990749e+01 -0.00000000e+00 -5.42232520e-04]

  day_night_flag :

  Day

  title :

  EMIT L2A Estimated Surface Reflectance 60 m V001

  granule_id :

  EMIT_L2A_RFL_001_20230401T203751_2309114_002

  Orthorectified :

  True

We can plot this information to see the spectra.

point.hvplot.line(x='wavelengths',y='reflectance',color='black').opts(title=f"Latitude: {point.latitude.values:.3f} Longitude: {point.longitude.values:.3f}")

3.4 Calculating CWC

As mentioned in the background we use the surface reflectance to estimate CWC. The unique spectral signatures allow identification and quantification based upon the wavelength-dependent absorption coefficients of liquid water. The EMIT mission has applied similar approaches to identify dust source minerals as well as methane point source emissions. The path length of liquid water absorption can be estimated by utilizing a least squares inversion to minimize the residuals between the EMIT reflectance and the Beer-Lambert model (Green et al.,2006), which relates the wavelength-dependent absorption to the path length the photons are traveling through the material. During the inversion, the path lengths are iteratively adjusted to match the modeled spectra to the EMIT reflectance within the water absorption feature region from 850 to 1100 nm.

First, define a Beer-Lambert Model function that returns the vector of residuals between measured and modeled surface reflectance.

# https://github.com/isofit/isofit/blob/main/isofit/inversion/inverse_simple.py#L514C1-L532C17
def beer_lambert_model(x, y, wl, alpha_lw):
    """Function, which computes the vector of residuals between measured and modeled surface reflectance optimizing
    for path length of surface liquid water based on the Beer-Lambert attenuation law.

    Args:
        x:        state vector (liquid water path length, intercept, slope)
        y:        measurement (surface reflectance spectrum)
        wl:       instrument wavelengths
        alpha_lw: wavelength dependent absorption coefficients of liquid water

    Returns:
        resid: residual between modeled and measured surface reflectance
    """

    attenuation = np.exp(-x[0] * 1e7 * alpha_lw)
    rho = (x[1] + x[2] * wl) * attenuation
    resid = rho - y

    return resid

We need some lab measurements of the complex refractive index of liquid water to obtain the wavelength-dependent absorption coefficients. They are calculated by taking four times the product of Pi and the imaginary part of the refractive index, divided by wavelength. The refractive index of liquid water per wavelength is provided by the k_liquid_water_ice.csv in the data folder. We can also preview this data to get a better understanding of the information we are using.

wp_fp = '../data/k_liquid_water_ice.csv'
k_wi = pd.read_csv(wp_fp)
k_wi.head()

	wvl_1	T = 22°C	wvl_2	T = -8°C	wvl_3	T = -25°C	wvl_4	T = -7°C	wvl_5	T = 25°C (H)	wvl_6	T = 20°C	wvl_7	T = 25°C (S)	Index
0	666.7	2.470000e-08	NaN	NaN	NaN	NaN	660.0	1.660000e-08	650.0	1.640000e-08	650.0	1.870000e-08	650.12971	1.674130e-08	0
1	667.6	2.480000e-08	NaN	NaN	NaN	NaN	670.0	1.890000e-08	675.0	2.230000e-08	651.0	1.890000e-08	654.63616	1.777420e-08	1
2	668.4	2.480000e-08	NaN	NaN	NaN	NaN	680.0	2.090000e-08	700.0	3.350000e-08	652.0	1.910000e-08	660.69347	1.939950e-08	2
3	669.3	2.520000e-08	NaN	NaN	NaN	NaN	690.0	2.400000e-08	725.0	9.150000e-08	653.0	1.940000e-08	665.27314	2.031380e-08	3
4	670.2	2.530000e-08	NaN	NaN	NaN	NaN	700.0	2.900000e-08	750.0	1.560000e-07	654.0	1.970000e-08	669.88461	2.097930e-08	4

fig, axs = plt.subplots(2,4, figsize=(15, 6),  sharex=True, sharey=True, constrained_layout=True)
axs = axs.ravel()
col_n = 0
for i in range(0, 7):
    x = k_wi.iloc[:, col_n+i]
    y = k_wi.iloc[:, col_n+i+1]
    axs[i].scatter(x, y)
    axs[i].set_title(y.name)
    col_n+=1
fig.supylabel('imaginary parts of refractive index')
fig.supxlabel('wavelength')
plt.show()

Now define a function to get the desired data from the csv file.

# https://github.com/isofit/isofit/blob/dev/isofit/core/common.py#L461C1-L488C26
def get_refractive_index(k_wi, a, b, col_wvl, col_k):
    """Convert refractive index table entries to numpy array.

    Args:
        k_wi:    variable
        a:       start line
        b:       end line
        col_wvl: wavelength column in pandas table
        col_k:   k column in pandas table

    Returns:
        wvl_arr: array of wavelengths
        k_arr:   array of imaginary parts of refractive index
    """

    wvl_ = []
    k_ = []

    for ii in range(a, b):
        wvl = k_wi.at[ii, col_wvl]
        k = k_wi.at[ii, col_k]
        wvl_.append(wvl)
        k_.append(k)

    wvl_arr = np.asarray(wvl_)
    k_arr = np.asarray(k_)

    return wvl_arr, k_arr

Lastly, to calculate CWC we define a function that uses least squares optimization to minimize the residuals of our Beer-Lambert Model and find a likely path length of liquid water.

# https://github.com/isofit/isofit/blob/main/isofit/inversion/inverse_simple.py#L443C1-L511C24
def invert_liquid_water(
    rfl_meas: np.array,
    wl: np.array,
    l_shoulder: float = 850,
    r_shoulder: float = 1100,
    lw_init: tuple = (0.02, 0.3, 0.0002),
    lw_bounds: tuple = ([0, 0.5], [0, 1.0], [-0.0004, 0.0004]),
    ewt_detection_limit: float = 0.5,
    return_abs_co: bool = False,
):
    """Given a reflectance estimate, fit a state vector including liquid water path length
    based on a simple Beer-Lambert surface model.

    Args:
        rfl_meas:            surface reflectance spectrum
        wl:                  instrument wavelengths, must be same size as rfl_meas
        l_shoulder:          wavelength of left absorption feature shoulder
        r_shoulder:          wavelength of right absorption feature shoulder
        lw_init:             initial guess for liquid water path length, intercept, and slope
        lw_bounds:           lower and upper bounds for liquid water path length, intercept, and slope
        ewt_detection_limit: upper detection limit for ewt
        return_abs_co:       if True, returns absorption coefficients of liquid water

    Returns:
        solution: estimated liquid water path length, intercept, and slope based on a given surface reflectance
    """
    
    # Ensure least squares is done with float64 datatype
    wl = np.float64(wl)
    
    # params needed for liquid water fitting
    lw_feature_left = np.argmin(abs(l_shoulder - wl))
    lw_feature_right = np.argmin(abs(r_shoulder - wl))
    wl_sel = wl[lw_feature_left : lw_feature_right + 1]

    # adjust upper detection limit for ewt if specified
    if ewt_detection_limit != 0.5:
        lw_bounds[0][1] = ewt_detection_limit

    # load imaginary part of liquid water refractive index and calculate wavelength dependent absorption coefficient
    # __file__ should live at isofit/isofit/inversion/
    
    
    data_dir_path = "../data/"
    path_k = os.path.join(data_dir_path,"k_liquid_water_ice.csv")
    
    #isofit_path = os.path.dirname(os.path.dirname(os.path.dirname(__file__)))
    #path_k = os.path.join(isofit_path, "data", "iop", "k_liquid_water_ice.xlsx")

    # k_wi = pd.read_excel(io=path_k, sheet_name="Sheet1", engine="openpyxl")
    # wl_water, k_water = get_refractive_index(
    #     k_wi=k_wi, a=0, b=982, col_wvl="wvl_6", col_k="T = 20°C"
    # )
    k_wi = pd.read_csv(path_k)
    wl_water, k_water = get_refractive_index(
        k_wi=k_wi, a=0, b=982, col_wvl="wvl_6", col_k="T = 20°C"
    )
    kw = np.interp(x=wl_sel, xp=wl_water, fp=k_water)
    abs_co_w = 4 * np.pi * kw / wl_sel

    rfl_meas_sel = rfl_meas[lw_feature_left : lw_feature_right + 1]

    x_opt = least_squares(
        fun=beer_lambert_model,
        x0=lw_init,
        jac="2-point",
        method="trf",
        bounds=(
            np.array([lw_bounds[ii][0] for ii in range(3)]),
            np.array([lw_bounds[ii][1] for ii in range(3)]),
        ),
        max_nfev=15,
        args=(rfl_meas_sel, wl_sel, abs_co_w),
    )

    solution = x_opt.x

    if return_abs_co:
        return solution, abs_co_w
    else:
        return solution

3.4.1 CWC of a Single Point

Now that we have all of the pieces, we can estimate the CWC of a single pixel using our invert_liquid_water function.

ewt = invert_liquid_water(point.reflectance.values,point.wavelengths.values)
print(f"EWT for ({point.longitude.values:.3f},{point.latitude.values:.3f}): {ewt[0]:.3f} cm")

EWT for (-120.353,34.540): 0.149 cm

3.4.2 CWC of a DataFrame of Points

We can also apply this function to the point data we selected in the previous notebook with the interactive plot.

Read in the csv file we created.

points_df = pd.read_csv("../data/emit_click_data.csv")
points_df

	id	x	y	381.00558	388.4092	395.81583	403.2254	410.638	418.0536	425.47214	...	2426.4402	2433.8303	2441.2183	2448.6064	2455.9944	2463.3816	2470.7678	2478.153	2485.5386	2492.9238
0	0	-120.569492	34.574552	0.018037	0.017798	0.017567	0.017348	0.017523	0.018065	0.018957	...	0.034822	0.034534	0.030492	0.031268	0.027023	0.026327	0.024637	0.023704	0.015793	0.013392
1	1	-120.568286	34.591054	0.012029	0.011957	0.011888	0.011826	0.012031	0.012486	0.013151	...	0.020698	0.022204	0.020115	0.019911	0.016761	0.017437	0.015805	0.016146	0.014353	0.012923
2	2	-120.371752	34.484274	0.020833	0.021323	0.021819	0.022317	0.023065	0.024095	0.025361	...	0.059581	0.058129	0.053620	0.053622	0.048783	0.044683	0.044064	0.041562	0.034080	0.028254
3	3	-120.181848	34.483303	0.024848	0.028178	0.034840	0.038940	0.043875	0.047232	0.050809	...	0.131800	0.135400	0.122960	0.118844	0.108201	0.112605	0.106271	0.101222	0.070507	0.059903
4	4	-120.407924	34.591054	0.040804	0.044085	0.048493	0.054733	0.059349	0.064003	0.067749	...	0.266081	0.263995	0.258724	0.252158	0.228920	0.226200	0.212124	0.181446	0.130708	0.111522
5	5	-120.438851	34.712298	0.015452	0.016034	0.016616	0.017204	0.017848	0.018621	0.019477	...	0.046427	0.046706	0.043317	0.044288	0.038928	0.037683	0.038662	0.035936	0.034616	0.033542
6	6	-120.512401	34.648230	0.044019	0.039991	0.039511	0.045176	0.045753	0.047011	0.048643	...	0.081432	0.080438	0.076009	0.076724	0.069166	0.069164	0.066233	0.063523	0.059582	0.056769
7	7	-120.350229	34.626389	0.018255	0.018679	0.019107	0.019550	0.020291	0.021369	0.022694	...	0.050578	0.047565	0.046092	0.046638	0.039017	0.036911	0.037903	0.033778	0.028663	0.028218
8	8	-120.355052	34.713269	0.017943	0.018290	0.018641	0.019003	0.019586	0.020435	0.021492	...	0.047776	0.048551	0.045117	0.045340	0.039024	0.038904	0.039077	0.035324	0.030073	0.029010
9	9	-120.243522	34.645318	0.021041	0.021242	0.021448	0.021661	0.022252	0.023214	0.024442	...	0.049286	0.047242	0.042362	0.045277	0.039692	0.040241	0.037825	0.034632	0.029134	0.029248

10 rows × 288 columns

Now create an array of wavelengths from the column names starting with the first wavelength value (column 3).

# Get wavelength values
wavelength_values = points_df.columns[3::].to_numpy()

Iterate by row through our dataframe, selecting the reflectance values and providing them to the invert_liquid_water function. Afterwards, add an CWC column to our dataframe.

# Creat empty list
ewt_values = []
# Iterate through rows
for _i in points_df.index.to_list():
    # Get reflectance array to pass to function
    rfl_values = points_df.iloc[_i,3::]
    # Use invert liquid water function and append results to list
    ewt_values.append(invert_liquid_water(rfl_values,wavelength_values)[0])    
# Add to our existing dataframe at Column Index 3
points_df.insert(3, "ewt", ewt_values)

points_df

	id	x	y	ewt	381.00558	388.4092	395.81583	403.2254	410.638	418.0536	...	2426.4402	2433.8303	2441.2183	2448.6064	2455.9944	2463.3816	2470.7678	2478.153	2485.5386	2492.9238
0	0	-120.569492	34.574552	2.773249e-01	0.018037	0.017798	0.017567	0.017348	0.017523	0.018065	...	0.034822	0.034534	0.030492	0.031268	0.027023	0.026327	0.024637	0.023704	0.015793	0.013392
1	1	-120.568286	34.591054	2.495512e-01	0.012029	0.011957	0.011888	0.011826	0.012031	0.012486	...	0.020698	0.022204	0.020115	0.019911	0.016761	0.017437	0.015805	0.016146	0.014353	0.012923
2	2	-120.371752	34.484274	1.639785e-01	0.020833	0.021323	0.021819	0.022317	0.023065	0.024095	...	0.059581	0.058129	0.053620	0.053622	0.048783	0.044683	0.044064	0.041562	0.034080	0.028254
3	3	-120.181848	34.483303	5.938623e-02	0.024848	0.028178	0.034840	0.038940	0.043875	0.047232	...	0.131800	0.135400	0.122960	0.118844	0.108201	0.112605	0.106271	0.101222	0.070507	0.059903
4	4	-120.407924	34.591054	4.001294e-11	0.040804	0.044085	0.048493	0.054733	0.059349	0.064003	...	0.266081	0.263995	0.258724	0.252158	0.228920	0.226200	0.212124	0.181446	0.130708	0.111522
5	5	-120.438851	34.712298	6.409762e-02	0.015452	0.016034	0.016616	0.017204	0.017848	0.018621	...	0.046427	0.046706	0.043317	0.044288	0.038928	0.037683	0.038662	0.035936	0.034616	0.033542
6	6	-120.512401	34.648230	2.607293e-10	0.044019	0.039991	0.039511	0.045176	0.045753	0.047011	...	0.081432	0.080438	0.076009	0.076724	0.069166	0.069164	0.066233	0.063523	0.059582	0.056769
7	7	-120.350229	34.626389	1.505599e-01	0.018255	0.018679	0.019107	0.019550	0.020291	0.021369	...	0.050578	0.047565	0.046092	0.046638	0.039017	0.036911	0.037903	0.033778	0.028663	0.028218
8	8	-120.355052	34.713269	1.379378e-01	0.017943	0.018290	0.018641	0.019003	0.019586	0.020435	...	0.047776	0.048551	0.045117	0.045340	0.039024	0.038904	0.039077	0.035324	0.030073	0.029010
9	9	-120.243522	34.645318	1.954748e-01	0.021041	0.021242	0.021448	0.021661	0.022252	0.023214	...	0.049286	0.047242	0.042362	0.045277	0.039692	0.040241	0.037825	0.034632	0.029134	0.029248

10 rows × 289 columns

For larger sets of point data we would want to do this in parallel. We can also calculate CWC for an area, where we definitly need to utilize parallel processing.

3.5 CWC Calculation of an ROI

In the previous notebook, we subset our region of interest and exported the file. Since the CWC calculation is computationally intensive, it can take a while to process large scenes, so its more efficient to do this spatial subsetting up front. We can use a function included in the ewt_calc.py module to calculate CWC on a cropped image, and create a cloud-optimized GeoTIFF (COG) file containing the results.

Set our input filepaths and output directory.

fp = '../../shared/2023-VITALS-Workshop-AGU/data/cropped/dangermond/rfl/EMIT_L2A_RFL_001_20230401T203751_2309114_002_dangermond.nc'

out_dir = '../data/'

roi_ds = xr.open_dataset('../../shared/2023-VITALS-Workshop-AGU/data/cropped/dangermond/rfl/EMIT_L2A_RFL_001_20230401T203751_2309114_002_dangermond.nc', decode_coords='all')
roi_ds

<xarray.Dataset>
Dimensions:           (wavelengths: 285, latitude: 244, longitude: 262)
Coordinates:
  * wavelengths       (wavelengths) float32 381.0 388.4 ... 2.486e+03 2.493e+03
    fwhm              (wavelengths) float32 ...
    good_wavelengths  (wavelengths) float32 ...
  * latitude          (latitude) float64 34.57 34.57 34.57 ... 34.44 34.44 34.44
  * longitude         (longitude) float64 -120.5 -120.5 -120.5 ... -120.4 -120.4
    elev              (latitude, longitude) float32 ...
    spatial_ref       int64 ...
Data variables:
    reflectance       (latitude, longitude, wavelengths) float32 ...
Attributes: (12/39)
    ncei_template_version:             NCEI_NetCDF_Swath_Template_v2.0
    summary:                           The Earth Surface Mineral Dust Source ...
    keywords:                          Imaging Spectroscopy, minerals, EMIT, ...
    Conventions:                       CF-1.63
    sensor:                            EMIT (Earth Surface Mineral Dust Sourc...
    instrument:                        EMIT
    ...                                ...
    spatialResolution:                 0.000542232520256367
    spatial_ref:                       GEOGCS["WGS 84",DATUM["WGS_1984",SPHER...
    geotransform:                      [-1.20992414e+02  5.42232520e-04 -0.00...
    day_night_flag:                    Day
    title:                             EMIT L2A Estimated Surface Reflectance...
    granule_id:                        EMIT_L2A_RFL_001_20230401T203751_23091...

xarray.Dataset

• Dimensions:
  • wavelengths: 285
  • latitude: 244
  • longitude: 262
• Coordinates: (7)
  • wavelengths
    (wavelengths)
    float32
    381.0 388.4 ... 2.486e+03 2.493e+03

    long_name :

    Wavelength Centers

    units :

    nm

    array([ 381.00558,  388.4092 ,  395.81583, ..., 2478.153  , 2485.5386 ,
           2492.9238 ], dtype=float32)

  • fwhm
    (wavelengths)
    float32
    ...

    long_name :

    Full Width at Half Max

    units :

    nm

    [285 values with dtype=float32]

  • good_wavelengths
    (wavelengths)
    float32
    ...

    long_name :

    Wavelengths where reflectance is useable: 1 = good data, 0 = bad data

    units :

    unitless

    [285 values with dtype=float32]

  • latitude
    (latitude)
    float64
    34.57 34.57 34.57 ... 34.44 34.44

    long_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    standard_name :

    latitude

    array([34.57432 , 34.573777, 34.573235, ..., 34.443642, 34.443099, 34.442557])

  • longitude
    (longitude)
    float64
    -120.5 -120.5 ... -120.4 -120.4

    long_name :

    longitude

    units :

    degrees_east

    axis :

    X

    standard_name :

    longitude

    array([-120.499254, -120.498711, -120.498169, ..., -120.358815, -120.358273,
           -120.357731])

  • elev
    (latitude, longitude)
    float32
    ...

    long_name :

    Surface Elevation

    units :

    m

    [63928 values with dtype=float32]

  • spatial_ref
    ()
    int64
    ...

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -120.49952471405297 0.0005422325202563338 0.0 34.57459066750936 0.0 -0.0005422325202563691

    [1 values with dtype=int64]

• Data variables: (1)
  • reflectance
    (latitude, longitude, wavelengths)
    float32
    ...

    long_name :

    Surface Reflectance

    units :

    unitless

    [18219480 values with dtype=float32]

• Indexes: (3)
  • wavelengths
    PandasIndex

    PandasIndex(Float64Index([ 381.0055847167969,  388.4092102050781,  395.8158264160156,
                  403.22540283203125, 410.63800048828125,  418.0535888671875,
                   425.4721374511719,  432.8927001953125,  440.3172607421875,
                   447.7427978515625,
                  ...
                   2426.440185546875,  2433.830322265625,   2441.21826171875,
                     2448.6064453125,  2455.994384765625,  2463.381591796875,
                   2470.767822265625,  2478.153076171875,   2485.53857421875,
                      2492.923828125],
                 dtype='float64', name='wavelengths', length=285))

  • latitude
    PandasIndex

    PandasIndex(Float64Index([ 34.57431955124923,  34.57377731872898,  34.57323508620872,
                  34.572692853688466, 34.572150621168205,  34.57160838864795,
                    34.5710661561277,  34.57052392360744, 34.569981691087186,
                  34.569439458566926,
                  ...
                   34.44743714150925,  34.44689490898899, 34.446352676468734,
                  34.445810443948474,  34.44526821142822,  34.44472597890796,
                   34.44418374638771,  34.44364151386745, 34.443099281347195,
                  34.442557048826934],
                 dtype='float64', name='latitude', length=244))

  • longitude
    PandasIndex

    PandasIndex(Float64Index([-120.49925359779284, -120.49871136527258, -120.49816913275234,
                  -120.49762690023208, -120.49708466771182, -120.49654243519157,
                  -120.49600020267131, -120.49545797015105, -120.49491573763079,
                  -120.49437350511054,
                  ...
                  -120.36261100268824, -120.36206877016798, -120.36152653764773,
                  -120.36098430512747, -120.36044207260721, -120.35989984008695,
                  -120.35935760756671, -120.35881537504645, -120.35827314252619,
                  -120.35773091000594],
                 dtype='float64', name='longitude', length=262))

• Attributes: (39)

  ncei_template_version :

  NCEI_NetCDF_Swath_Template_v2.0

  summary :

  The Earth Surface Mineral Dust Source Investigation (EMIT) is an Earth Ventures-Instrument (EVI-4) Mission that maps the surface mineralogy of arid dust source regions via imaging spectroscopy in the visible and short-wave infrared (VSWIR). Installed on the International Space Station (ISS), the EMIT instrument is a Dyson imaging spectrometer that uses contiguous spectroscopic measurements from 410 to 2450 nm to resolve absoprtion features of iron oxides, clays, sulfates, carbonates, and other dust-forming minerals. During its one-year mission, EMIT will observe the sunlit Earth's dust source regions that occur within +/-52° latitude and produce maps of the source regions that can be used to improve forecasts of the role of mineral dust in the radiative forcing (warming or cooling) of the atmosphere.\n\nThis file contains L2A estimated surface reflectances and geolocation data. Reflectance estimates are created using an Optimal Estimation technique - see ATBD for details. Reflectance values are reported as fractions (relative to 1). Geolocation data (latitude, longitude, height) and a lookup table to project the data are also included.

  keywords :

  Imaging Spectroscopy, minerals, EMIT, dust, radiative forcing

  Conventions :

  CF-1.63

  sensor :

  EMIT (Earth Surface Mineral Dust Source Investigation)

  instrument :

  EMIT

  platform :

  ISS

  institution :

  NASA Jet Propulsion Laboratory/California Institute of Technology

  license :

  https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/

  naming_authority :

  LPDAAC

  date_created :

  2023-04-06T02:26:24Z

  keywords_vocabulary :

  NASA Global Change Master Directory (GCMD) Science Keywords

  stdname_vocabulary :

  NetCDF Climate and Forecast (CF) Metadata Convention

  creator_name :

  Jet Propulsion Laboratory/California Institute of Technology

  creator_url :

  https://earth.jpl.nasa.gov/emit/

  project :

  Earth Surface Mineral Dust Source Investigation

  project_url :

  https://emit.jpl.nasa.gov/

  publisher_name :

  NASA LPDAAC

  publisher_url :

  https://lpdaac.usgs.gov

  publisher_email :

  lpdaac@usgs.gov

  identifier_product_doi_authority :

  https://doi.org

  flight_line :

  emit20230401t203751_o09114_s000

  time_coverage_start :

  2023-04-01T20:37:51+0000

  time_coverage_end :

  2023-04-01T20:38:03+0000

  software_build_version :

  010610

  software_delivery_version :

  010610

  product_version :

  V001

  history :

  PGE Run Command: {python /beegfs/store/emit/ops/repos/emit-sds-l2a/spectrum_quality.py /tmp/emit/ops/emit20230401t203751_emit.L2AReflectance_20230405t233102/output/emit20230401t203751_rfl /tmp/emit/ops/emit20230401t203751_emit.L2AReflectance_20230405t233102/output/emit20230401t203751_rfl_quality.txt}, PGE Input Files: {radiance_file=/beegfs/store/emit/ops/data/acquisitions/20230401/emit20230401t203751/l1b/emit20230401t203751_o09114_s000_l1b_rdn_b0106_v01.img, pixel_locations_file=/beegfs/store/emit/ops/data/acquisitions/20230401/emit20230401t203751/l1b/emit20230401t203751_o09114_s000_l1b_loc_b0106_v01.img, observation_parameters_file=/beegfs/store/emit/ops/data/acquisitions/20230401/emit20230401t203751/l1b/emit20230401t203751_o09114_s000_l1b_obs_b0106_v01.img, surface_model_config=/beegfs/store/emit/ops/repos/emit-sds-l2a/surface/surface_20221020.json}

  crosstrack_orientation :

  as seen on ground

  easternmost_longitude :

  -119.71545648976226

  northernmost_latitude :

  34.7990749308955

  westernmost_longitude :

  -120.992414074966

  southernmost_latitude :

  33.77100207248943

  spatialResolution :

  0.000542232520256367

  spatial_ref :

  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

  geotransform :

  [-1.20992414e+02 5.42232520e-04 -0.00000000e+00 3.47990749e+01 -0.00000000e+00 -5.42232520e-04]

  day_night_flag :

  Day

  title :

  EMIT L2A Estimated Surface Reflectance 60 m V001

  granule_id :

  EMIT_L2A_RFL_001_20230401T203751_2309114_002

roi_ds.sel(wavelengths=850,method='nearest').hvplot.image(x='longitude',y='latitude',cmap='viridis',geo=True, tiles='ESRI', frame_width=720,frame_height=405, alpha=0.7, fontscale=2).opts(
    title=f"Dangermond ROI - RFL at 850 nm", xlabel='Longitude',ylabel='Latitude')

Use the calc_ewt function to calculate CWC of the cropped image. This function will also create a COG file containing the CWC results. We can also specify the number of CPUs to use manually with a n_cpu argument, or leave it blank to use all but one of the available CPUs. If we set the return_cwc argument to true, the function will also return the COG.

This will take some time, about 5 minutes, because we’re doing the calculation for roughly 63,000 pixels.

%%time
ds_cwc = calc_ewt(fp, out_dir, return_cwc=True)

2023-12-07 15:11:32,965 WARNING services.py:1996 -- WARNING: The object store is using /tmp instead of /dev/shm because /dev/shm has only 67108864 bytes available. This will harm performance! You may be able to free up space by deleting files in /dev/shm. If you are inside a Docker container, you can increase /dev/shm size by passing '--shm-size=3.05gb' to 'docker run' (or add it to the run_options list in a Ray cluster config). Make sure to set this to more than 30% of available RAM.
2023-12-07 15:11:34,128 INFO worker.py:1664 -- Started a local Ray instance. View the dashboard at 127.0.0.1:8265 

CPU times: user 688 ms, sys: 315 ms, total: 1 s
Wall time: 5min

ds_cwc

<xarray.Dataset>
Dimensions:      (latitude: 244, longitude: 262)
Coordinates:
  * latitude     (latitude) float64 34.57 34.57 34.57 ... 34.44 34.44 34.44
  * longitude    (longitude) float64 -120.5 -120.5 -120.5 ... -120.4 -120.4
    elev         (latitude, longitude) float32 ...
    spatial_ref  int64 ...
Data variables:
    cwc          (latitude, longitude) float64 nan nan nan nan ... nan nan nan
Attributes: (12/13)
    flight_line:            emit20230401t203751_o09114_s000
    time_coverage_start:    2023-04-01T20:37:51+0000
    time_coverage_end:      2023-04-01T20:38:03+0000
    easternmost_longitude:  -119.71545648976226
    northernmost_latitude:  34.7990749308955
    westernmost_longitude:  -120.992414074966
    ...                     ...
    spatialResolution:      0.000542232520256367
    spatial_ref:            GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84...
    geotransform:           [-1.20992414e+02  5.42232520e-04 -0.00000000e+00 ...
    day_night_flag:         Day
    title:                  EMIT Estimated Equivalent Water Thickness (EWT) /...
    granule_id:             EMIT_L2A_RFL_001_20230401T203751_2309114_002

xarray.Dataset

• Dimensions:
  • latitude: 244
  • longitude: 262
• Coordinates: (4)
  • latitude
    (latitude)
    float64
    34.57 34.57 34.57 ... 34.44 34.44

    long_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    standard_name :

    latitude

    array([34.57432 , 34.573777, 34.573235, ..., 34.443642, 34.443099, 34.442557])

  • longitude
    (longitude)
    float64
    -120.5 -120.5 ... -120.4 -120.4

    long_name :

    longitude

    units :

    degrees_east

    axis :

    X

    standard_name :

    longitude

    array([-120.499254, -120.498711, -120.498169, ..., -120.358815, -120.358273,
           -120.357731])

  • elev
    (latitude, longitude)
    float32
    ...

    long_name :

    Surface Elevation

    units :

    m

    [63928 values with dtype=float32]

  • spatial_ref
    ()
    int64
    ...

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -120.49952471405297 0.0005422325202563338 0.0 34.57459066750936 0.0 -0.0005422325202563691

    [1 values with dtype=int64]

• Data variables: (1)
  • cwc
    (latitude, longitude)
    float64
    nan nan nan nan ... nan nan nan nan

    long_name :

    Canopy Water Content

    units :

    g/cm^2

    _FillValue :

    -9999

    array([[nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           ...,
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan],
           [nan, nan, nan, ..., nan, nan, nan]])

• Indexes: (2)
  • latitude
    PandasIndex

    PandasIndex(Float64Index([ 34.57431955124923,  34.57377731872898,  34.57323508620872,
                  34.572692853688466, 34.572150621168205,  34.57160838864795,
                    34.5710661561277,  34.57052392360744, 34.569981691087186,
                  34.569439458566926,
                  ...
                   34.44743714150925,  34.44689490898899, 34.446352676468734,
                  34.445810443948474,  34.44526821142822,  34.44472597890796,
                   34.44418374638771,  34.44364151386745, 34.443099281347195,
                  34.442557048826934],
                 dtype='float64', name='latitude', length=244))

  • longitude
    PandasIndex

    PandasIndex(Float64Index([-120.49925359779284, -120.49871136527258, -120.49816913275234,
                  -120.49762690023208, -120.49708466771182, -120.49654243519157,
                  -120.49600020267131, -120.49545797015105, -120.49491573763079,
                  -120.49437350511054,
                  ...
                  -120.36261100268824, -120.36206877016798, -120.36152653764773,
                  -120.36098430512747, -120.36044207260721, -120.35989984008695,
                  -120.35935760756671, -120.35881537504645, -120.35827314252619,
                  -120.35773091000594],
                 dtype='float64', name='longitude', length=262))

• Attributes: (13)

  flight_line :

  emit20230401t203751_o09114_s000

  time_coverage_start :

  2023-04-01T20:37:51+0000

  time_coverage_end :

  2023-04-01T20:38:03+0000

  easternmost_longitude :

  -119.71545648976226

  northernmost_latitude :

  34.7990749308955

  westernmost_longitude :

  -120.992414074966

  southernmost_latitude :

  33.77100207248943

  spatialResolution :

  0.000542232520256367

  spatial_ref :

  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

  geotransform :

  [-1.20992414e+02 5.42232520e-04 -0.00000000e+00 3.47990749e+01 -0.00000000e+00 -5.42232520e-04]

  day_night_flag :

  Day

  title :

  EMIT Estimated Equivalent Water Thickness (EWT) / Canopy Water Content (CWC)

  granule_id :

  EMIT_L2A_RFL_001_20230401T203751_2309114_002

Plot CWC of our ROI.

ds_cwc.hvplot.image(x='longitude',y='latitude',cmap='viridis',geo=True, tiles='ESRI', frame_width=720,frame_height=405, alpha=0.7, fontscale=2).opts(
    title=f"{ds_cwc.cwc.long_name} ({ds_cwc.cwc.units})", xlabel='Longitude',ylabel='Latitude')

Contact Info:

Email: LPDAAC@usgs.gov
Voice: +1-866-573-3222
Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹
Website: https://lpdaac.usgs.gov/
Date last modified: 11-28-2023

¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.